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ABSTRACT 

We perform axisymmctric relativistic magnetohydrodynamic (MHD) simulations to investigate the 
acceleration and coUimation of jets and outflows from disks around compact objects. Newtonian 
gravity is added to the relativistic treatment in order to establish the physical boundary condition 
of an underlying accretion disk in centrifugal and pressure equilibrium. The fiducial disk surface 
(respectively a slow disk wind) is prescribed as boundary condition for the outflow. We apply this 
technique for the flrst time in the context of relativistic jets. The strength of this approach is that it 
allows us to run a parameter study in order to investigate how the accretion disk conditions govern the 
outflow formation. Substantial effort has been made to implement a current-free, numerical outflow 
boundary condition in order to avoid artiflcial coUimation present in the standard outflow conditions. 
Our simulations using the PLUTO code run for 500 inner disk rotations and on a physical grid size of 
100x200 inner disk radii. The simulations evolve from an initial state in hydrostatic equilibrium and an 
initially force-free magnetic fleld configuration. Two options for the initial field geometries are applied 
- a hourglass-shaped potential magnetic field and a split monopolc field. Most of our parameter runs 
evolve into a steady state solution which can be further analyzed concerning the physical mechanism at 
work. In general, we obtain coUimated beams of mildly relativistic speed with Lorentz factors up to 6 
and mass-weighted half-opening angles of 3-7 degrees. The split-monopole initial setup usually results 
in less coUimated outflows. The light surface of the outflow magnetosphere tends to align vertically - 
implying three relativistically distinct regimes in the flow - an inner sub-relativistic domain close to 
the jet axis, a (rather narrow) relativistic jet and a surrounding subrelativistic outflow launched from 
the outer disk surface - similar to the spine-sheath structure currently discussed for asymptotic jet 
propagation and stability. The outer subrelativistic disk-wind is a promising candidate for the X-ray 
absorption winds that are observed in many radio-quiet AGN. The hot winds under investigation 
acquire only low Lorentz factors due to the rather high plasma-/3 we have applied in order to provide 
an initial force-balance in the disk-corona. When we increase the outflow Poynting flux by injecting 
an additional disk toroidal field into the outfiow, the jet velocities achieved are higher. These flow 
gains super-magnetosonic speed and remains Poynting flux dominated. 

Subject headings: accretion, accretion disks - ISM: jets and outflows - MHD - galaxies: active — 
galaxies: jets - relativity 



1. INTRODUCTION 

Astrophysical jets emanate from sources spanning a 
huge range in energy output or length scale - among them 
young stellar objects (YSO), stellar mass compact ob- 
jects as X-ray binaries or /i-quasars, or the powerhouses 
of some active galactic nuclei (AGN) which host a super- 
massive black hole. In particular for radio- loud quasars, 
for which synchrotron emission dominates the radio spec- 
trum, relativistic jets are a generic feature. Due to the 
omnipresent angular momentum conservation, mass ac- 
cretion to all of these objects features a disk structure 
around the central mass. It is commonly believed that 
jets are launched as disk winds, which are fur ther accel- 
erated and coUimated by ma gnetic forces (see 'Blandford 
fc Payne| (ll982t ; jPudritz fc N orman (1983); Camcnzind 



1986U;|He^||ll)l)7| ; |Hey^aerts fc Norman| ( j2003p ;|Pu- 
dritz et al. ( |2007| . Relativistic jets may gam further en- 



ergy by interaction with the black hole magnetosphere 
( [Blandford fc ZnajeklpTT) [Ghosh fc Abramowicz||l997 
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|Komissarov"2005) . 

The magnctohydrodynamic (MHD) self-coUimation of 
non-relativistic jets has been proven in general by time 
dependent sim ulations ( Ustyugova et al ||l995 Ouyed fc 



Pudritz 19971 and have been investigated in further de 
tail considerin g addition al physic al effe cts as magnetic 
diffusivity by J endt fc Cemeljic"] (20021, a variation in 
Ouyed fc Pudritz ( 1999 )7^on-axisymmet ric instabilities 



m the launching region ( )Ouyed et al. 2003[) , or a variation 
in the mass flow profile or the magnetic tie 



leld geometries 

([Fendt 2006 Pudritz et al.[[2006t , or the influence of a 



central magnetic held ( i^'endt|2009 Matsakos et al.|2008" 



In the case of relativistic jets the efficiency ot MHD 
self-coUimation is under debate. The main reason is the 
existence of electric fields which are negligible for non- 
relativistic MHD and which are commonly thought to 
have a net de-coUim ating effect on the jet. Essentially, 



Chiueh et al. (19911 have demonstrated the current car- 
rying relativistic jet can be highly coUimated. However, 
the actual structure of these jets still remains unclear - 
mainly due to the need for simplifying assumptions to 
solve the corresponding set of MHD equations. 
So far, a variety of theoretical models have been devel- 
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oped for the case of s elf- similar ]eis ( Li et al.|1992 
topoulos|[T994| [T995I IVlahakis fc K6mgl||!^()0i]n ^ 
et al.||2006p, aith ough it seems clear that relativity 



Con- 



eliani 

does 

not obey self-similarity. Fully 2.5D theoretical solutions 
for the internal magnetic jet str ucture could be obtained 



by neglec ting matter inertia ( ^Fendt 1997a Fendt & Mem- 
ola 2001). These force- free solutions for the held struc- 



ture can in principle be coup led to the dynamical wind 



solution along the field lines (Fendt fc Camenzind 



Fendt & Greiner 2001 l-'endt h Ouyed'SUD^ 



\ 1997a[ ) obtamed solutions tor the internal jet force bal- 



1996 



'tendt 



ance m Kerr metric with an asymptotically cylindrical 
jet emerging from a disk-like structure around the cen- 
tral rotating black hole. The shape of the coUimating 
jet boundary was obtained as result of the internal force 
equilibrium, in particular considering the regularity con- 
dition along the jet outer light surface. 

Timc-dcpcndcnt simulations of relativistic MHD jet 
formation have been performed considering a general rel- 
ativistic metric, including also the evolution of the un- 
derlying accretion disk. Early - seminal - simulations 



did last for a few inner disk rotations only (Koide et al. 
fl998, 1999), which is sufficient time to demonstrate the 
launching of an outflow, but hardly sufficient in order 
to investigate the long-term dynamical evolution of the 
emerging jet. 

More recent simulations were able to follow several 100 
disk rotations and show the formation of a so-called fun- 



zon and the inner disk radius (jDe Villiers et al. 


2005 


McKinney & Narayan||2007| 


i'chekhovskoy et al. 


2008 


McKinney & Blandford| 2009 


1. These simulations indi- 



degree of coUimation. However, the funnel flow achieves 
Lorentz factors of up to 50. General relativistic MHD 
simulations are also able to determine the interrelation 
betwe en jet for mation and the Blandford-Znajek mech a- 



nism (McKinney 2005: Komissarov & McKinney 2007) 



Ultra-relativistic MHD simulations of accel erating and 
coUimating jets have been presented by Komis sarov et ah] 
([2007), spanning over a huge range of length scale and 
providing jets of large Lorentz factor F ~ 10. Their sim- 
ulations, however, did not start from the very base of 
the jet - the accretion disk, but at some fiducial bound- 
ary above the equatorial plane. Since the jet has been 
launched already with super-escape speed, gravity has 
not been considered. The jet flow has been confined 
within a rigid wall of predefined shape which naturally 
affects the opening angle of the MHD jet nozzle and thus 
jet collimation and acceleration. 

The focus of our present paper is i) to concentrate 
on the formation and acceleration of a relativistic MHD 
jet right from the launching area the accretion disk sur- 
face, ii) to investigate the (self-) collimation of relativistic 
MHD jets under the infiuence of de-collimating electric 
forces and an "open" boundary condition for the outfiow 
iii) to consider gravity as en essential gradient to provide 
a realistic disk boundary condition in equilibrium, iv) to 
run long-term simulations lasting more than 1000 inner 
disk rotations until the jet reaches steady state, v) to 
concentrate on MHD disk jets as disks are the natural 
origin for the mass load for AGN jets. 

The outline of this paper is as follows: In section [2] 
we discuss the concepts of ideal special relativistic mag- 



netohydrodynamics in the perspective of jet formation. 
Section [3] is devoted to the initial- and boundary con- 
ditions of the numerical simulations, whose results are 
shown in section ID We conclude in section [5] 

2. CONCEPTS OF RELATIVISTIC MHD JETS 

It is well k nown that relativistic jets must b e strongly 
magnetized ( |Michel|p69| |Camenzind]|1986b| p| |l993| ). 
This simply reflects the fact that the lower tne mass 
flux, the more electro-magnetic energy (Poynting flux) 
can be transferred into high kinetic energy per unit mass. 
Relativistic MHD is also limited for very strong magne- 
tization as then the MHD assumption can be violated 
since a sufficiently large amount of electric charges is 
lacking which are needed to drive the electric current sys- 
tem. Such a situation might arrise in the ultra-relativistic 
regime of pulsar-winds but is unlikely for the disk- winds 
investigated here. 

This paper deals with the time-dependent formation 
of relativistic MHD jets by using the special relativistic 
MHD modul e of the PLUTO code provided by [Mignonej 
et al. (2007) and applying a Newtonian description of 



gravity. 

2.L Relativistic MHD equations 

The relativistic MHD module of PLUTO solves the 
system of special relativistic conservation laws. In a co- 
variant formulation, the equations follow naturally as a 
set of hyperbolic equations. It is solved for energy and 
momentum conservation 



(1) 



a^T"'' = 

of an ideal magnetized fiuid 

= {ph + b^)u"u^ +(^P+ ^b^^ g'^f^ - h'^b^ (2) 
with the specific plasma enthalpy 

/i=^- + l, (3) 
7 - 1 p 

the isotropic gas pressure p, density p (both in the local 
rest-frame), four-velocity (u") = (F,F/3)-^, velocity (3 — 
v/c, Lorentz factor F = (1 — /3^)~°-^ and the magnetic 
field pseudo vector 

b''^~\e''^^'upF,s. (4) 

Assuming infinite conductivity, thus vanishing electric 
fields in the rest-frame of the plasma F^-f^up = 0, the 
homogenous Maxwell equation 



d„ ^ 



(5) 



can be written solely in terms of the magnetic four- vector 

= ^e"'^''*F^5 = b^uf^ - b'^u"' (6) 

and for the field vector components it follows^ 

B' = = b'u° ~ b°u' (7) 
E'^e'^^Vu^. (8) 



^ following the convention that Greek indices run from to 4 
whereas Latin indices go from 1 to 3 
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Equation [8] represents the ideal MHD condition E = 
—(3 X B and is the reason why all electric fields can be 
eliminated from the equations. The magnetic four-vector 
turns out as 6° = B'u'; b' = (B' + b°u')/u°. The con- 
servation of the Faraday tensor given by Eq. [6] results in 
the non-rclativistic (ideal) induction equation and the 
solenoidal condition V • B — 0. Mass conservation is 
guaranteed by the continuity equation 



(9) 



We apply a polytropic equation of state for the gas with 
the polytropic index 7 — 5/3. 

2.2. Gravity in special relativity 

The outcome of MHD simulations is mainly deter- 
mined by its boundary conditions and therefore requires 
great care in describing the proper physical state of inter- 
est. Since in our simulations the jet is considered to be 
launched as a wind from a rotating disk, it is essential to 
take into account a proper disk model as boundary condi- 
tion. For the disk boundary we choose a (sub-) Keplerian 
rotatio n profi le and a hydrostatic pressure gradient (see 
section 



3.1.11 



This choice implies a hot, "puffed up" 
corona with a sound speed ~ 2/3 w^. This is the main 
reason why we have added gravity to our special rela- 
tivistic treatment. The direct impact of gravity on the 
jet dynamics is marginal as the outflow is accelerated to 
super escape speed quickly. 

A fully self-consistent relativistic treatment of gravity 
would imply a general relativistic approach. However, we 
are not interested in the region very close to the horizon 
of the central object, but in the long-term dynamics and 
evolution of a disk wind into a relativistic jet. We can 
therefore neglect general relativistic effects in our simu- 
lation domain. 
We apply a softened gravitational potential 



GM 
R + rs 



(10) 



with a softening length of = 1/3 that may be related 
to the to the Schwarzschild radius of a non-rotating black 
hole'^. The corresponding acceleration reads 



a = -V(f) = -- 



GM 



(11) 



{R. + rs)^R 

and hence instead of solving equation [T] we solve 

d^T""^ = (12) 

with the four force density {f^) — Tp{a ■ v, a)-^ as a local 
source term on the right-hand side. This is incorporated 
in PLUTO as a "body force" (Fa) using the infrastruc- 
ture of the code. 

Omission of softening would lead to numerical errors 
(due to the unresolved steep gradients in the potential 
close to the origin), piling up to produce artificial ac- 
celeration along the spine of the jet close to the axis. 
Softening is clearly a compromise avoiding the singular- 
ity (by limiting the required resolution) on little cost of 



3 The capital _R : 
throughout this work. 



\/r^~+~z^ denotes the spherical radius 



realism. Ano ther choice could be the wel l-known pseudo- 
potential by Paczynsky & Wiita ( 1980 1 which has just 
the negative softening 0pn/ = —GM/{R — rs). For the 
cylindrical geometry of our choice, the singularity would 
become even more problematic, complicating the setup 
a great deal. 



2.3. Relations in axisymmetric MHD 

The region of jet formation may be fairly well approx- 
imated in axisymmetry. In fact, non-axisymmetric dis- 
tortions may actually hinder the formation of powerful 
jets as probably demonstrated by the existence of a va- 
riety of strongly magnetized, rapidly rotating accretion 
disk systems which, however, do not exhibit jets (e.g. 
cataclysmic variables or most pulsars). 

Under the assumed symmetry in a cylindrical coordi- 
nate system, the magnetic field vector can be written as 



B = B, 



(13) 



where B^ can now be an arbitrary function of r and 
z, as the solenoidal condition translates to V • Bp = 0. 
The stream function ^{r, z) = (l/27r) /dS • Bp = rA^ 
measures the magnetic fiux through the surface area S 
and follows from the toroidal component of the vector 
potential. 



= V X Aa = V X = -V^- X ew, (14) 

r r 



For the electric field, 
E = B X /? gives 



the ideal MHD condition 



-BpU 



BpU 



(15) 



in terms of the so called angular velocity of the field line 
— {v^ — VpB^/Bp) /r, or the so-called hght cylinder 
radius of a field line tl = c/H.^ . The direction of the 
electric field is given by n = Bp /Bp x and is perpen- 
dicular to the magnetic fiux surface ^'(r, z). The poloidal 
Poynting fiux S = (c/47r)E x B^ simplifies to 



S=- 



An ri^ 



An 



(16) 



2.3.1. Perpendicular and parallel force-balance 

The processes leading to fiow coUimation can be iden- 
tified directly from the (steady-state) trans-fi eld force- 
balance equation (Chiueh et al. 1991 Appl fc Camei>] 
zind 1993). Here, we adopt the notation of the latter 
paper when investigating the coUimation behavior in the 
quasi steady-state time domain of our simulations. The 
curvature K = n (Bp • V) Hp/ Bp of a fiux surface ^'(r, z) 

* This is the radius at which the hypothetical angular velocity 
of a field line supercedes the speed of light 
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results from the summation of perpendicular forces, 




(17) 



47rc2 



-rpVj 



where we have added the coUimating compo- 
nent of the gravit ationa l force. For the ease 

we label the terms as 
1 the 



of use in section |4.2.1 



inch; F^cii - ^ol 7 -^grav) 



; ^pbp ; -^pbphi ; -^p ■} -^piii _ 

order of their appearance in EqUT^ The (poloidal) 
Alfven Mach number M is relativistically defined as 



El ■ 



(18) 



The gradient Vj^ = n • V is projected perpendicular to 
the magnetic flux surfaces , and thus along the (inward 
pointing) electric field. The light surface of a magne- 
tosphere is located where ri^{^) — ri^{r,z) = c/^l^{^) 
^, it hence depends on the flux-geometry as well as on 
the rotation profile fi^. Each fiux surface / magnetic 
field line crosses the light surfa ce at most once (see also 
the discussion in Fendt (1997b)) Some field lines ^{r,z) 
never cross the light surface, indicating an asymptotic 
radius roo(^') < J'l,*- For these field lines relativistic 
effects due to rotation (electric fields) are less impor- 
tant. For others, the asymptotic radius is roo(^') > tl,*. 
The light surface constitutes a critical point of the sta- 
tionary axisymmetric wind equation only in the mass- 
less limit in which it is identical to the modified polo idal 
Alfven surface M\ = 1 — (r^/rL)^ ( Camenzind|l986a|bj ) . 
However, it is essential to note that at the light surface 
the dynamical behavior of the poloidal magnetic pressure 
term changes - the force changes sign. This leads to the 
existence of three dynamically different regimes in the 
asymptotic (coUimated) region of a relativistic jet (see 
Fig.fl]). In region I, for all field lines r^oi^) < ^'l,* corre- 
sponding to (1 — (r/rL)^) > 0, and, thus, a de-coUimating 
magnetic pressure term. Field lines in region II do cross 
their light cylinder, and, since r > tl, the magnetic pres- 
sure term acts as coUimating for r > tl. Field lines 
in Region III never reach their light cylinder, and here 
the magnetic pressure term is de-coUimating again. The 
slope of the outer part of the light surface critically de- 
pends on the dynamics and the magnetic field structure 
of the outfiow in the very inner part. 
Similarly, the forces due to the electric field E = 



r/ri^Bp (second last term in Eq. 17 1 scale with the rela- 
tive position to the light surface, hence they are impor- 
tant in regi on I I of the jet formation region only. 

Equation [17] together with Fig.fl] once more demon- 
strates the need to resolve the whole acceleration and 
coUimation region of a jet in radial and vertical direction. 

^ Thus, the Ught surface consists of the points of intersection 
between the field lines with their corresponding light cylinder 



Only when the light surface is taken into account self- 
consistently, the proper force-balance is applied along 
and across the flow. 




Fig. 1. — The different dynamical regimes in special relativistic 
disk-winds. Region 1 and 111 stay sub-relativistic, while regi on II 
is relativistic, i.e. electric forces are not negligible (see Sect. |2.3~T] 
for a discussion). 



Similarly one can derive the parallel-field force equa- 
tion, it becomes: 



An 



w / Br 




(19) 



-rpV||c 



with the necessary definitions V|| = Bp/_BpV and 
K\\Bp = Bp/Bp (Bp • V) Bp. We see how the change in 
the Mach- number is mediated by the interplay of tension- 
, pressure-, pinch-, centrifugal- and gravitational- accel- 
eration. Electric fields (pointing in the perpendicular 
direction) don't contribute and the equation reduces to 
the Newtonian case: In steady-state, there is no electric 
acceleration! 

A number of self-similar approaches to the relativis 
tic jet formation have beed pubhshed (eg 
fc Konigl] (l2003|) 



Vlahakis 

While the self-similar ansatz is a 
powerfull and highly successfull tool to solve the non- 
relativistic MHD problem (starting with the Blandford- 
Payne solution), we believe that using self-similarity for 
relativistic MHD jets is problematic. 

We note that neither the light surface nor the rela- 
tivistic Alfven surface obeys a self-similar structure. It 
is well known that forcing self-similarity into the rela- 
tivistic MHD equations constrains the rotation law for 
the magnetosph ere ^ l^(r) oc r~^. (see also discussion in 
Li et al.] ( |1992[ ); p] (]l993|). This is a major difference 
to the non-relativlslic self-similar approach. We further 
note that also the scaling for the electric field depends 
on the radial position of the light surface (see Eq. 15 1. 



This is, however, of uttermost importance for the struc- 
ture of relativistic magnetospheres as the electric field 
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forces pla y a leading role in the trans-field force-balance 
(equation 17 1. Similar arguments hold for the inner light 



integral part in models of the X-ray features of AGN 



surfaces around Kerr black holes or the geometry of the 
black hole ergosphere. 

We therefore believe that a steady-state self-similar rel- 
ativistic MHD approach is intrisically inconsistent with 
the relativistic characteristic of the flow. 

2.3.2. field line constants 

Stationary axisymmetric MHD flows conserve the fol- 
lowing five quantities along the magnetic fiux- function ^P. 
From the iso-rotation law together with the ideal MHD 
condition follows the rest mass energy-flux per magnetic 
induction. 

pUp 



k = /c(«') = 
and the iso-rotation parameter 



1 



Br, 



v„ — - 



(20) 



(21) 



(often interpreted as angular velocity of the field lines). 
In absence of shocks the (pseudo-) entropy 

Q = ^ = Qi^) (22) 
is conserved as well as the angular momentum flux 



27rfcc 



ru^ = /(*) 



(23) 



and the flux ratio of total energy to rest-mass energy, 



S + IC + M+T + g 
M 



m(*) 



(24) 



where we identify the individual terms as (purely) ki- 
netic energy flux K. = (T ~ l)p Up, rest-mass energy flux 
Ai = p Up, thermal energy flux T = F^j^rP ^pi ^nd 
gravitational energy flux Q = p (p Up, respectively. The 



cold, asymptotic limit of (24 1 is particularly of interest, 
it reads 

p = Tia+l) (25) 

where a = S/{1C + M) is the customarily defined mag- 
netization parameter - the ratio of Poynting to kinetic 
fiux. This simple relation provides a theoretical max- 
imum for the Lorentz-factor T* = p, when the entire 
electromagnetic energy is converted into kinetic energy. 

The essential point in the quest for relativistic jets is 
to find a highly energetic disk solution with values of p 
beyond the anticipated Lorentz-facto r. Previous stud- 
i es obt aining highly relativistic jets by |Komissarov et aL] 
( 2007 1 do not start from a realistic disk solution but in- 



ject the jet material with an artificial rotation profile, a 
high injection speed Vp > 0.5c, and a density low enough 
to obtain a high energy flux 12 < /imax < 18. In the 
present study we are aiming to improve on the problems 
just mentioned by applying a physical boundary condi- 
tion as a Keplerian disk-corona in equilibrium. 

2.4. Accretion disk coronae 

It is our ambition to connect the wind solutions to 
the ambiance of a realistic accretion disk. In our sim- 
ulations, the flow originates in the high entropy atmo- 
sphere called a corona. Optically thin coronae are an 



(e.g. [M ushotzky et al. 1993) and /i-Quasars (e.g. Nowak 
et al .p002| |Markotf et al.||2003| ). 



while Compton cooling can provide the observed spec- 
tra, the heating mechanism is not easily found. Just as 
in the case of the sun, the coronal heat cannot directly 
be transfered from the colder photosphere/accretion disk 
(according to the second law of thermodynamics) and the 
nature of vertical energy transport is an active field of re- 
search. External irradiation of fiared disks by the central 
object (or cen tral dis k) is certainly present in a multitude 
of objects (see Czerny et al.|2008[ for a review) but might 
not be the primary energy source. Among the most 
promising mechanisms we should hi ghlight magnetic re- 
connec tion heating as proposed by Haardt fc Maraschi 

pgr). 

Between the mid-plane and the coronal point of injec- 
tion in our simulations, ideal MHD can not provide a 
realistic picture. In order for an accretion disk to work, 
a torque of viscous or magnetic origin has to be exerted 
onto the material. Additionally, the fiux-freezing con- 
straint of ideal MHD must be relaxed since it would lead 
to an accumulating magnetic pressure that ultimately 
stops the accretion process. Studies modelling both the 
accretion motion and the super Alfvenic jet based on a 
stationary self-similar approach (e.g . Wardle & Koenigl 



([1993 
globa 



Ferreira fc Pelletier| ( |1993| ; ^ ( |1995] )) rely on 
scale magnetic fields and ad-hoc assumptions on 



the viscosity and (anisotropic) magnetic diffusion. In or- 
der for these models to be stable against strong magnetic 
compression on the o ne side and the magneto r otational 
instab ility (MRI; se e 'Balbu s fc Hawley] ( |1998| ) on the 
other, Ferreira & Pellcticr ( 1995| require equipartition 
for the the rmal and magnetic pressure. |Casse fc Fer- ' 



reira 



( 2000 1 demonstrated the importance of heating for 

the niass-loading or the jet. Time-dependend numer- 
ical simulations of these magnetized accretion ejection 



structures (MAES) were presente d by Casse fc Keppens 



( |2002[ |2004[ ); [Zanni et ah] ([2007') and adopt a fixed in 
time anomalous resistivity profile in order to connect the 
two dynamical regimes. 

Although the stationarity of the aforementioned simu- 
lations is possibly hampered by numerical diffusion and 
low resolution effects, MAES with global-scale fields are 
to date the most successful models to create coUimated 
outfiows. 

It is now widely believed that the source of viscosity 
and resistivity in weakly magnetized disks is the turbu- 
lence seeded b y the MRI. Local, str atified shearing box 
simulations by Miller fc Stone ( 2000 1 suggest a quenching 
of the MRI in the strongly magnetized coronal region, as 
the magnetic scale height exceeds the thermal one. This 
is in contrast to the equipartition fields proposed for the 
MAES. Magnetic buoyancy of the large-scale fields even- 
tually created by a turbulent dynamo could then pro- 
vide the coronal heating. Typically, the MRI results in 
toroidally dominated coronal fields with i3| > lOB^ and 
opening the field lines towards a topology favorable for 
wind acceleration remains a challenge. In the context 



of the solar co rona this process is discussed by Wang fc 



Sheeley[([2003[). 
Hy restraining parts o f the accretion disk structure. 



von Rekowski et al. ( 2003 1 could achieve outfiows in time- 
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dependend simulations of disk-corona structures where 
the magnetic field is sustained by a mean-field dynamo. 
Analytic models of this turbulent disk-corona-outflow 
connection are however in still its in fancy (see the dis- 
cussion by Kuncic fc Bicknell| ( 2004 1 and attempts by 
packman fc Fcssah^ ( ,2009J ). 

3. MODEL SETUP FOR THE MHD SIMULATIONS 

With the aforementioned considerations we choose the 
following model for our investigation. A global-scale 
poloidal field favorable of wind acceleration is adopted. 
Whether it is advected by the accretion flow or created by 
an underlying dynamo is not of our concern. The jet base 
resembles a corona in the sense that it is hot (electron 
temperature lO^K), has no mechanism of cooling, is 
non-turbulent (no viscosity) and highly ionized (inflnite 
conductivity). We choose a Keplerian rotation profile 
for the field lines. The fiow starts with sub-escape veloc- 
ity and we investigate sub-magnetosonic injection where 
mass loading is determined by the internal dynamics as 
well as mass-fluxes imposed by the boundary condition. 
We perform axisymmctric special relativistic MHD sim- 
ulations of jet formation for a set of different magnetic 
field geometries and field strengths. In the following we 
discuss the numerical realization of our problem. 

3.1. Boundary conditions 

Given the 2.5 dimensional nature of the problem, three 
geometrical boundaries have to be prescribed. These are 
the inlet boundary along z = from which material 
is injected into the domain (inflow) and the two outer 
boundaries at r = rend and z = Zcnd where we expect 
material to leave the computational domain (outflow). 
The boundary condition along r — (Rbeg) follows from 
cylindrical symmetry. Figure [2] gives an overview of the 
different regions. 



3.1.1. Injection boundary (Zteg) 

Pursuing the aim to follow the acceleration of an disk 
wind from as close to the accretion disk as possible, 
we start with a sub slow-magnetosonic wind. We are 
hence free to choose four constraining boun dary condi- 
tions w ithout overdetermining the system (see Bogovalov 
(19971 and Appendix [B] for more details) 



Our choice is to fix the toroidal electric field component 
= 0. This suppresses the evolution of the bound- 
ing poloidal magnetic field and is realized by requiring 
Vp||Bp. The field line iso-rotation is kept constant in 
time and follows a Keplerian rotation law fi^ cx r~^-^. 
Unless specified otherwise, the boundary condition starts 
out initially in a force-free state with zero toroidal field 
f2^ = v^{r)/r corresponding to a disk in hydrodynamic 
equilibrium. This is an essential ingredient as - within 
stationary ideal MHD - ^l^ just equals the mid-plane an- 
gular velocity of the material uj{r). Relaxation of the in- 
finite conductivity constraint would, however, lead to an 
inequality f2^(r) < uj{r) owing to the diffusion of mag- 
netic field. For these reasons il^ should closely follow the 
expected disk rotational profile and should be limited by 
the maximal velocity in the mid-plane, typically at the 
inner edge of the disk located at r = 1 . 

A radial force-equilibrium along the whole boundary is 
enforced by balancing the centrifugal and pressure sup- 
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Fig. 2. — Sketch of the different regimes of our grid and boundary. 
In both directions we set 20 equidistant cells in [0, 1]. Then follows 
a stretched grid until we add five equidistant cells one unit-radius 
before the outflow boundary. For r £ [0, 1], the hydrostatic corona 
is fixed to minimize the influence of the central region on the disk- 
wind. 

port against gravity via the sub-Keplerianity of the ro- 
tation ^Jx = '^4>iT = ^)I'Vk that also determines the in- 
let density, where vk is the circular velocity that alone 
sustains against gravity at r = 1. A more convenient 
parametrization is in terms of the relative temperature 



^M = (7-l)(l 



x)/x- We investigate two cases 
2 /3 and a version with e = 1/6 



- a hot corona with e 
(X = 0.5 and x = 0-8) 

If we interpret r = 1 as the innermost stable circular 
orbit (ISCO) around a black hole, uk is a measure of 
the black hole spin. In the case of a Schwarzschild black 
hole it is — 0.6c while we choose the scaling velocity 
v,^(r = 1) = 0.5c for convenience. 

The inner disk-edge is numerically difficult to model 
because of the transition to the infiow of the disk-wind 
and the steep gradients in gravity. Within r < 1, the 
so-called plunge region, a physical solution would allow 
for (radial and vertical) accretion onto the central ob- 
ject. Since the dynamics in this area would then require 
a general relativistic treatment which we cannot provide 
in this context, we simply minimize the dynamical effect 
of this region by freezing the hydrostatic solution initially 
present in the domain. Other auth ors h ave assumed a 
thin fu nnel flow along the axis (e.g. [Krasn opolsky ct aL 
(|1999 or a dded an internal sink-cell (lik e |Casse fc Kep- 
pens(2002l) to circumvent this problem 'i'he transition 
between the "inner corona" and disk-wind is smoothed 
via the Fermi step function F(x) = (1 -I- e'y^~^')l'^-'^)~'^ . 
Rotational support is thus turned on by the setting 



Pikr,z) = - 



1 



1 - xF^{R) 
v4r) = ^VKF{R)R-°-'. 

Density given by equation [26] and the coronal pressure p 



-iR + rs)'/^'-''\ (26) 
(27) 
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constitute the third and fourth fixed in time conditions. 

We emphasize that is is not possible to specify both 
injection velocity and density-profile and thus the mass- 
fiux for sub-(magneto)sonic fiows, as this is determined 
by the sonic point. Therefore we match the vertical ve- 
locity to the domain via d^Vz = 0, while the radial 
component follows from the Erj, = condition. We limit 
the injection speed by the local slow magnetosonic speed 
in the case when the velocity just above the boundary 
becomes trans-sonic. This provides the fifth constraint 
needed in that case. 

With the induction of a toroidal magnetic field com- 
ponent in the jet, the rotatio nal velocity needs to be ad- 
justed in order to satisfy Eg. [21] 



W0 — rfl^ 



Br, 



(28) 



as we apply the condition cJz-B^ — dzBr = (while 
then follows from V • B = 0). 

By letting the jet solution alone determine Poynting- 
and mass- fiux, we loose control over the energy fiux 
parameter /x and the limiting asymptotic Lorentz factor 
r*. It will rather be a consequence of the MHD under 
the constraints we have given, while we have used our 
freedom to provide a boundary most closely resembling 
a realistic hot disk corona. A graphical summary of the 
disk wind boundary conditions is shown in Fig. [3] 

In order to extend the parameter space towards 
higher /i, we also investigate cases where we have over- 
determined the boundary conditions by specifying the 
mass-flux through 



similar to 



e.g. 



Vz{r, 0) = Wi„jW0(r, 0) 
Ouyed & Pudritz] p97| ; 



(29) 



Fendt & 



Cemeljic (2002). Another parameter run adopts a 1/r 
profile tor the toroidal magnetic field B^{r) = —riF{r)/r 
in order to specify the Poynting flux with an additional 
parameter rj. These runs and the influence of the ov er- 
determination is discussed separately in section |4.3.2[ 
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Fig. 3. — Profiles of the fixed in time variables for the inlet in hy- 
drodynamic equilibrium. Here we give constraints on Vl^ , p,p, 
{E^ = not shown). Parameters are ijk = 0.5, e = 2/3. The 
thin dotted line is the Fermi step function used to smooth those 
variables experiencing a sharp transition at the inner disk radius 
r = 1. 



3.1.2. Outflow boundaries (Rend, Zgnd) 



The standard outflow boundary conditions for many 
numerical codes are zero-gradient conditions, which are 
usually sufficient as the plasma velocities are constrained 
to be outward-pointing. 

In the case of sub fast-magnetosonic outflows, this 
strategy is unfortunately insufficient as the flow inside of 
the domain will depend on the flow beyond the bound- 
ary via the incoming characteristics. Just as for the inlet 
boundary, the now missing information has to be sup- 
plied by constraints that describe best the physical con- 
ditions downstream of the boundary. In the case of an 
outflow, the conditions leading to an untampered flow are 
however impossible to know a priori. A way to circum- 
vent this unphysical feedback is to avoid any causal con- 
tact by moving the boundary far away such that the char- 
acteristics will not enter the domain of interest within the 
simulated time. 

When considering a boundary outside of causal con- 
tact "very far away" , we estimate for Alfven waves to 
travel over lO'^ scale radii within the anticipated simula- 
tion time. The computational effort of such huge grids 
does not allow a large parameter study at the current 
time and we must leave this option for future endeav- 
ours. 

In the absence of a substantially better solution, zero- 
gradients are used for the primitive variables except for 
magnetic flelds for which this simple approach leads to 
artificial electric currents implying an inward-pointed 
Lorentz-force. Especially for low plasma-/? this may re- 
sult in a devastating artificial coUimation - preventing 
any steady state to establish and artificially coUimating 
t he outflow increasin gly thin with time. 

([1999| have performed a system- 



Ustyugova et al. 



atic study comparing ditierent approaches for outflow 
conditions including a (toroidal) force-free condition 
jpl |Bp = and a more sophisticated version including an 
additional numerical factor that needs to be determined 
a posteriori. For the outflow conditions in our simula- 
tions we instead recover the magnetic field components 
by imposing constraints on the poloidal (jr — —dzB^, 
jz — r~^drrB^) and toroidal (j^ = dzBr — drBz) electric 
currents. For the toroidal magnetic field (poloidal elec- 
tric current) we radially extrapolate the expected 1/r 
law of a marginal jz at the radial end (Rend)- Using 
dzB^ = allows to specify jr at the upper end of the 
domain (Zend)- Concerning the poloidal magnetic field 
components we implement a current-free boundary con- 
dition by enforcing = 0. This is a novel approach 
designed to minimize spurious effects of coUimation. We 
convinced oursevels that boundary effects have only a 
marginal effect on the solution by varying the grid-size 
and geometry. For a detailed discussion and comparison 
of various outflow conditions we refer to appendix \K\ 

We note that, as a further complication, a fully rela- 
tivistic version for a force-free or force-balance boundary 
conditions would also need to take into account electric 
forces. We have estimated the impact of such an upgrade 
and found that due to the geometry of our outflow (in 
particular the location of the light surface) it would play 
a minor role and is thus not worth the effort to imple- 
ment. 
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3.2. Initial conditions 

As initial state we prescribe a force-free coronal mag- 
netic field, F°-^jp = 0, together with a gas distribu- 
tion in hydrostatic equilibrium. Both is essential in or- 
der to avoid artificial relaxation processes caused by a 
non-equilibrium initial condition. We apply a polytropic 
equation of state p = Kp^ with a "classical" polytropic 
index of 7 = 5/3 since our flows are always cold when 
compared to the rest-mass. To further strengthen this 
choice , we p erformed a comparison simulation with the 



Taubl ( 1948 1 equation of state as described by Mignone 
et al.| ( [2UU5| which produced an identical jet once the hot 



shock has passed through. The constant K is determined 
by the radial force-balance of the inlet. 

For the initial magnetic fleld conflguration we apply 
two different geometries. Field conflguration A is a po- 
tential field of hour glass shape as appl i ed by |Ouyed fc| 



Pudritz 



(|l997l and iFendt & Cemeljicl (|2002| with the 



magnetic field components 



Br 



1 



1 



{z + Zd) 



(r2 + (z + zd)2)'/' 
1 

(r2 + (z + ^d)2)'/'' 
corresponding to a vector potential 
1 



v/r-2 + {z + ZdY - (z + Zd) 



(30) 



(31) 



(32) 



in cylindrical coordinates with i?,. = —dzA^ and B^ = 
r^^dr rA^. The dimensionless disk thickness with 
(zd + -z) > is introduced to avoid kinks in the field 
distribution for z < (the ghost zones) and we choose 
Zd = 1 for convenience. 

Our other option for the initial m agnetic field (c onfig- 
uration B) is the "split monopole" (Sakurai 1987| with 
the magnetic field components 



Br — 



B, 



(r2 + (z + z,)2)3/2 

Z + Zd 



(r2 + (z + z<j)2) 



3/2 ■ 



(33) 



(34) 



In the split- monopole setup the parameter Zd defines the 
offset of the fiducial center of the monopole from the grid 
origin, and adjusts the initial angle of the field lines with 
respect to the disk-surface. We either adopt an angle of 
6 = 85° or 9 = 77° for the field line passing through 
r — \. Similar to Eq. [32j the split monopole field can be 
described by a vector potential 



Aa, = 



1 



1 - 



Z + Zd 



(r2 + (z + Zd)2) 



2^l/2 



(35) 



The fields are scaled to satisfy the the choice of the 
plasma-/? 



/3: 



(36) 



r=l,z=0 



at the inner disk radius. 

It should be kept in mind that plasma-/? largely varies 
along the disk boundary. In configuration A, the profile 



/3(r) monotonically decreases until for large radii it is 
/3(r) cx r^^-^ leading to a magnetically dominated outer 
corona. In the split-monopole, /3(r) decreases first to a 
minimum value (at r*{9 = 77°) « 5 and r*{9 = 85°) « 
15) and increases for large radii according to I3{r) cx r^'^ 
leading to thermal dominance. 

In summary, for our injection boundary condition we 
are left with the following five dynamical parameters. 



(uK,/3,e,Winj,'7), 



(37) 



where strictly speaking we are only allowed to choose the 
first three when launching sub-slow. An overview of the 
simulations performed in this parameterization is shown 
in Tab. [21 

3.3. Numerical grid and physical scaling 

We use a numerical grid of 512 x 1024 cells applying 
cylindrical coordinates. Onward from the inner region 
(r < l,z < 1), which is resolved with 20 x 20 equidis- 
tant cells, we apply a stretched grid with the element 
size increasing by a factor of < 1.005. This leads to a 
domain size of (r x z) = (100 x 200)ri corresponding to 
(300 X 600) if n = 3rs (see sketch in figure |2|. Stag- 
gered magnetic fields tr eated via constrained transport 
( [Balsara fc Spicer|[l999 1 are used to ensure V • B = 0. 

Because of the constraints imposed on the cell aspect 
ratio by the zero-current boundary (appendix A|, we set 
the last five grid cells to be equally spaced witnmaximal 
aspect ratios < 3/1. 

The dimensionless nature of our simulations allows for 
various astrophysical interpretations. We provide a phys- 
ical scaling of simulation variables (marked with a prime) 
in the following paragraph. 

Since velocities are given in terms of the speed of light 
(c' = 1), relativistic simulations are in need of only two 
additional scales. The simulation variables are connected 
to their physical counterparts via 

v = v'c; 1 = 1'Iq; t ^ t' ~ t' If) j vq; p — p' Po (38) 

P = p'po=pW; B = B'So = BV4^Poc2. (39) 

If we assume a Schwarzschild black hole as central body, 
we may set the spatial scale = Sr^, equating the inner 
disk radius with the ISCO. Then it becomes 



Wo = 3 X 10^°cm s"^ 



M. 



/o = 9 X lO^cm 



in = 3 X lO^^s 



(40) 
(41) 

(42) 



Assuming a physical outfiow mass-loss rate in terms of 
the Eddington limited accretion rate M = O.OlMedd we 
can provide a scale for the density by comparison to the 
mass loss rate of the simulation M' 



po = 6 X 10 



, 1 



M, 

Mo 



g cm 



(43) 



where we applied a radiative efficiency of 7y* = 0.1. The 
scaling of pressure and magnetic fields then follows as 
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TABLE 1 
Fiducial scaling 



M, 


l/V 


t/t' 


P/P' 


viv' 


B/B' 




[cm] 


\A 


[gem 


[g s— cm--*-] 


[Gauss] 


10« 


9 X lO^^ 


3 X 10^ 


1.8 X IQ-l'^ 


1.5 X 10^ 


1.4 X 10^ 


10 


9 X 10'^ 


3 X 10-"' 


1.8 X 10-9 


1.5 X 10^2 


4.4 X 10^ 



Note. — Scaling for simulation WA05 {M' — 32.67) assuming a 
physical mass loss rate of M — l%Medd with an efficiency of r?* — 0.1. 



P0=5 X 10 



14 1 (M, 



M' \Mq 



— 1 —9 

cm ^ s 



(44) 



Bo = 8 X 10' M 



7 f 



-0.5 



Gauss. (45) 



Under these considerations, the only remaining scahng 
parameter is the mass of the compact object M,. Ne- 
glecting additional physical processes as radiation pres- 
sure or radiative cooling leaves us with a scale-free model 
that can be applied to any disk- wind launched jet around 
compact objects. Table [T] provides a fiducial scaling for 
a microquasar with M, = lOM© and for an AGN with 
M, = 10* M0. The scale- free nature becomes obvious if 
we recall the res t-frame teni peratures for an ideal gas, 



T 



p_ 

p'k-B 



{^i)inp (Anile 



19891, with the mean molecular 



weight (fi), the proton mass mp, and the Boltzmann con- 
stant k^. For ionized hydrogen (/i) = 0.5 one would find 
temperatures of T = {p' / p') 5.45 x lO^^K, while for an 
electron-positron plasma {fj, ~ 1/2000) the temperatures 
are lower by three orders of magnitude. These ultra- 
high scaling temperatures are only reached in the very 
inner corona between the central object and the inner 
disk radius. As we do not intend to follow the dynamics 
here, this does not really pose a problem. In principle, 
once p'/p' ^ 1, the equation of state transcends towards 
7 = 4/3 according to a Synge-gas. In the jet, thermal 
pressure quickly looses importance and the temperatures 
are significantly lower. 

4. RESULTS AND DISCUSSION 

We now present the results of our numerical simula- 
tions considering the formation of relativistic MHD jets 
from accretion disks. Each simulation consumed approx- 
imately 48 hours on 16 processors. The overall goal is 
to test whether the paradigm of MHD self-coUimation of 
non-relativisti c jets es t ablish e d from numerical simula - 
tions fUstyugovaet al.| (|1995|; |Ouyed fc Pudrit'il(|1997| ; 



Krasnopolsky et al. ( |l999| ; |Fendt fc Cemeljic| ( |2002| also 
holds in the relativistic case. 

4.1. Overall evolution of the outflow 

The initial evolution of the disk corona is governed by 
the propagation of toroidal Alfven waves launched due to 
the rotation of the field line foot-points. The initial force- 
free magnetic field structure is adapted to a new dynamic 
equilibrium according to a rotating wind magnetosphere. 

A wind is launched from the disk boundary and is con- 
tinuously accelerated driving a shock front through the 
initial hydrostatic corona and sweeping this material out 
of the computational domain (Fig. HI). The disk wind 



evolves into a coUimated outflow of super-magnetosonic 
speed. Along the symmetry axis the hydrostatic initial 
condition is very well preserved. Once the bow shock has 
passed through the domain, the jet mass flux declines to a 
value which is solely governed by the internal outflow dy- 
namics and the injection boundary conditions. Similarly, 
the post-shock magnetic field distribution follows as well 
from the internal outflow dynamics and has in principle 
little in common with the initial setup. Certain combina- 
tions of boundary conditions for mass flux and magnetic 
fleld will result in a quasi-stationary^ state of the outflow 
evolution (see next section). From this point onwards we 
can start our investigations of coUimation and accelera- 
tion. In this paper we concentrate on analysis when the 
flow has reached a quasi-steady state. We usually ter- 
minate our simulations after 500 inner disk-rotations P, 
while a quasi-steady state is established over most of the 
domain after about 200 rotations. 

Figure |4] shows the time evolution for two exemplary 
simulations with an initial hourglass-shaped potential 
field distribution (case A) and a split-monopole field 
distribution (case B), each for the parameter choice 
(/3,«;^,e) = (1,0.5,2/3). 

The figure shows the Lorentz factor, the poloidal mag- 
netic field lines, poloidal electric current fiow lines and 
the critical MHD surfaces. In addition the light surface 
is drawn. 

Phenomenologically, the solutions form a magnetic 
nozzle with, depending on the disk fiux distribution, con- 
siderable difference in the width, but comparable final 
opening angles of the fast component. A broader ini- 
tial field distribution (case B) also results in broader 
and faster winds where the material originating from 
the inner disk is more effectively thinned out. In anal- 
ogy to hydrodynamic nozzles, the flow reaches the slow- 
magnetosonic speed directly above the throat. CoUima- 
tion happens mainly before the fast-magnetosonic sur- 
face is reached. Afterwards, the opening angle of a given 
fleld line is approximately conserved. 

Of particular interest is the electric current distribu- 
tion (shown for the time step T/P — 250). The electric 
current distribution is a consequence of the dynamical 
evolution of the outflow and therefore a direct outcome 
of the disk boundary magnetic flux profile and the Kep- 
lerian field line rotation. 

In general, the electric current leaves the outer disk 
to return within the fast component of the outflow. It 
is expected to enter the inner disk and then flow radi- 
ally outwards closing with the outgoing current. Such 
butterfly-shaped circuits are expected in Keplerian disks 
while the ir pl ays a leading role in the disk-jet feedback 
( Ferreira]|1997 1. A positive radial electric current in the 



disk-corona supports accretion by braking the disk ma- 
terial due to its magnetic torque jV x similar to a 
Barlow- wheel^. 

^ We denote the dynamieal state as quasi stationary as due to 
Keplerian disk rotation the disk outflow a large radii has evolved 
for a considerably lower number of disk rotations. Therefore, a 
slight change in 1;he dynamical state of the outer outflow can be 
expected after another few outer, thus 1000s of inner rotations. 

^ However, this region is not resolved within our numerical 
domain, as it is located below our injection boun dary as part of the 
underl y ing non-ideal MHD accretion disk. See |Casse fc Keppensj 
poMj l; |Zanni et a£\ pOOTf for non-relativistic simulations ot the 
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TABLE 2 

Parameter summary of our disk- wind simulations. 



ID 


Top 


/9 


e 


Vinj 


V 


Remarks 


Tmax 


/^max 




"^p , max 


Tjet 


M 


WAOl 


A 


0.2 


2/3 


var. 


var. 




1.23 


1.33 


18.84 


0.54 


20.26 


23.77 


WA02 


A 


1 


2/3 


var. 


var. 




1.27 


1.37 


11.47 


0.58 


21.86 


26.62 


WA03 


A 


2 


2/3 


var. 


var. 




1.26 


1.34 


10.69 


0.57 


22.21 


24.15 


WBOl 


B 


1 


2/3 


var. 


var. 


e = 77° 


1.33 


1.41 


8.27 


0.64 


24.19 


16.48 


WB02 


B 


1 


2/3 


var. 


var. 


e = 85° 


1.29 


1.33 


8.19 


0.60 


21.27 


15.66 


WA04 


A 


0.2 


1/6 


var. 


var. 




1.27 


1.38 


13.40 


0.58 


21.14 


36.08 


WA05 


A 


1 


1/6 


var. 


var. 




1.25 


1.33 


9.84 


0.57 


22.77 


32.67 


WA06 


A 


2 


1/6 


var. 


var. 




1.25 


1.32 


9.34 


0.57 


22.92 


29.27 



Note. — Columns arc from left to right: simulation ID; initial magnetic field distribution (Top): potential field (A) or split monopole (B); 
plasma-/?; accretion disk temperature parameter e; injection speed Vi^j as a fraction of Vfp{r)) toroidal disk-field scaling ri; specific remarks; to the 
right we show values of the steady state, the maximum Lorentz factor Fmax collimation degree ^, maximal poloidal velocity Vp,max, jet radius rjet 
and the total mass flux out of the domain M. 



The inchnation between the poloidal current vector 
and the magnetic field line indicates the direction of (de- 
) coUiniating magnetic forces acting on the flow. When 
the inclination becomes less than 90°, the Lorentz force 
jp X changes from collimation to de-coUimation. This 
can be clearly seen in the snapshots at T/P = 250 of the 
case A simulation where actual field lines are indicated in 
white and initial field lines in red. In the actual field dis- 
tribution, the field lines are somewhat pushed away from 
the surface jp _L Bp (this is also where magnetic acceler- 
ation is most effective). For case B this happens beyond 
the light surface. As a result, both electric and magnetic 
forces deflect the flow towards the disk boundary which 
leads to a highly unstable layer just above the outer disk. 
We will provide an in-depth analysis including all forces 
acting on the flow in section |4.2.1[ 

The locations of the characteristic surfaces are signa- 
tures for the MHD flow. Depending on the initial mag- 
netic flux di stribution (c ase A,B) and the mass flux pro- 
flle (see also Fendt (2006)), this location may vary a great 
deal. In our case B simulations, we generally observe sur- 
faces which leave the domain in radial direction (parallel 
to the disk surface) . For the case A simulations these sur- 
faces tend to "coUimate" leaving the domain in vertical 
direction. The latter implies a two-layered structure of 
the jet - a central super-fast magnetosonic jet surrounded 
by a sub-Alfvenic outflow. This is an interesting aspect 
for observational modeling and for stability analysis of 



faces at a certain height (here at z = Zm), 



2007 



rma 



2005; Mizuno et al. 



sheath-spine jets (fP ushkarev et al. 

IKovalev et ai:|2007( |Hardee||20o7nBeskin fe Mokh 
2009| . The broad wind launched from the outer 



regions ot the disk has much lower velocities, decreasing 
continously with increasing launching-radius. For exam- 
ple, the terminal velocity of the flow originating from 
rfp > 32 of the case A simulations drops below 0.2c, con- 
sistent with the X-ray abso rption features observed in a 
mounting number of AGN ( Cappi|2006 Turner & Miller 
20091. In principle, our dynamical models can provide 
basic ingredients (e.g. flow geometries and velocity gra- 
dients) for the mod eling of spectral hne p rofiles of disk 
winds dKnigge et al . 1995; Sim et al.||2008 l. 

Following F^nHt J 2006| ) , we may define an average col- 
limation degree f of the outflow measured as the fraction 
of vertical and radial mass flux through equal-area sur- 

disk-jet interaction 



./2 



rmTVrp\r„ dz 



(46) 



(47) 



Similarly, we deflne a mass flux weighted jet radius, 

_ Jo"" rTv^p\^^ dr 
/J^"* Tw^pI^^ dr 

The corresponding values for ^ and rjot derived for 
Zm — 200 at the upper end of the domain and for time 
t/P = 500 are given in Tab.|2] along with the maxi- 
mum Lorentz factor Fmaxi the maximum poloidal veloc- 
ity Vp,max, and the total mass flux M. Figure [H] shows 
the time evolution of these quantities in the top panel. 
In general, we observe that the collimation degree ^ is 
the most sensitive tracer for secular trends among the 
observables mentioned. In the lower panel, we show the 
evolution of jet-power in the individual energy channels 
leaving the computational domain (radial and vertical). 
After the re-configuration of the initial stationary-state 
to the dynamical solution, the partitioning of energies is 
completed at around t/P = 100. Thermal energy-flux 
peaks when the hot bow-shock passes through the upper 
boundary. Far away from the central object, gravita- 
tional and thermal energy-flux are negligible The inte- 
grated energy-flux is dominated by rest-mass, reflecting 
the fact that only the inner component reaches signiflcant 
Lorentz-factors. The balance between Poynting and ki- 
netic flux is of particular interest. Figure [5] shows merely 
the end result of the spatial conversion history with the 
remaining electromagnetic energy S above the purely ki- 
netic part /C. More detailed insight into how this is es- 
tablished is provided in the following section using an 
individual field line. 

4.2. Stationary state analysis 

Simulations starting from an initial field distribution A 
evolve into a quasi-stationary flow solution after about 
200 inner disk rotations. Figure [7] shows our reference 
simulation WA04 at time t/P — 250, including an en- 
larged subgrid of the innermost area of the domain. 

Steady state solutions are helpful to understand the 
flow structure for a number of reasons. Firstly, by using 
MH D cons ervation laws, the conserved quantities (see 
Sect. 2.3.2) allow to identify the momentum and energy 
channels of the flow during acceleration and collimation. 
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Fig. 4. — Formation of a relativistic MHD jet. Shown is the Lorentz factor (color gradient) in terms of logio(r ~ 1) at the time of 25, 
250 and 500 (left to right) inner disk rotations. Top: Field distribution case A, hourglass-shape potential field (run WA02). Bottom: Split 
monopole initial field distribution case B (run WBOl). Shown are poloidal magnetic field lines (solid white lines); the critical MHD surfaces 
(solid black lines); the light surface r ■ Q = c (solid black line). For time T/P = 250 electric current flow lines are added (solid green lino). 
Time step T/P = 500 also show the initial field lines for comparison (dashed white lines). 



Secondly, by using the force-balance equations |17|19| we 
may identify the leading forces on the material along the 
outflow. Thirdly, the cross-check for conserved quanti- 
ties provides another test for the quahty of our setup 
and the numerical approach. A secondary indicator of 
stationarity is the alignment of poloidal velocities with 
the poloidal magnetic field lines, = 0. Figure 17] shows 
corresponding velocity vectors confirming this picture. 

This is confirmed by checking in detail the complete 
set of integrals of motion of the MHD-fiow k, fi^, Q,l, ^ 



as defined in equations [20]- 24 Figure [6] shows the rela- 



tive deviation of these quantities from their average value 
along a given field line after t/P — 500. The integrals 
are conserved within 1%-accuracy already right above 
the injection boundary - clearly demonstrating the qual- 
ity the choice of our numerical setup, in particular the 
injection boundary conditions carefully constructed from 
an equilibrium of Keplerian rotation and gas pressure. 

Due to the differential rotation law, the number of Ke- 
plerian rotations t/P{r) scales with radius as t/P{r) = 
{t/ P)r~^/'^ , implying that at the end of our simulations 
{t/P — 500), we have performed roughly one rotation 
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Fig. 5. — Time evolution of characteristic quantities. Top panel: 
After an initial adjustment till t/P ^ 200, mass flux M, jet ra- 
dius rjet, coUimation degree ^, maximal Lorentz-factor Fmax and 
poloidal velocity I'p.max cease to evolve. Lower panel: Power in 
the individual energy chanels out of the domain. Thermal power 
(T) peaks when the shock reaches the upper boundary and is neg- 
ligible otherwise. The total outgoing power (labeled accordingly) 
is dominated by rest- mass (A/) and Poynting flux (5). Also shown 
are the gravitational (G) and (purely) kinetic (K) contributions 
(simulation run WA02). 

at r = 64 and half a rotation at r = 100. Nonetheless 
the integrals of motion for the field line rip — 64 are 
conserved within 0.1%. 
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Fig. 6. — Field line constants (conserved quantities). Shown is 
the deviation from the average value along the fleld line rooted at 
rfp = 2 (simulation run WA02). 
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Fig. 7. — Logarithmic (rest-frame) density of the stationary flow 
(simulation run WA04). Shown are poloidal magnetic fleld lines 
(solid white), electric current flow lines (solid black), characteristic 
MHD surfaces (various dot-dashed green) , surface of escape veloc- 
ity (dotted green), light surface (solid green). Arrows in the top 
plot indicate the velocity fleld. The bottom flgure is an enlarged 
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4.2.1. Collimating and accelerating forces 

In this section we identify the forces responsible for 
jet acceleration and colhmation applying the st eady-s tate 
parallel and transversal force-equilibrium Eqs. 



Fig.|8] compares these forces for a number of refer- 
ence simulations (WBOl, WA02, WA05) along a field 
line rooted at rfp = 2. As check for consistency, we also 
show the gradient of the Mach number a = Bp/AnW ^^A'P 
which just coincides with the summation of the parallel 
forces, indicating a steady state (see yellow solid and 
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black dashed line). 

In general, the outflow starts with sonic speed and is 
first launched by thermal pressure in the hot disk corona, 
respectively the centrifugal force in the colder version. 
Until the Alfven point, the Lorentz- force of the poloidal 
electric current (-F'pbphi + ^pinch) is the main magnetic 
driver. Ultimately the poloidal tension (-Fcurv) keeps the 
acceleration up even above the fast surface. 

Concerning the transverse force, we reproduce the ex- 
pected sign-change of the curvature (tension) force (first 
coUimating until the Alfven surface, de-coUimating be- 
yond) and the poloidal pressure force (de-coUimating un- 
til the light cylinder, coUimating beyond). For the cross- 
field balance, we observe the following three regimes: 

Just on top of the inlet, the main de-coUimating forces 
besides poloidal magnetic pressure are thermal pressure 
in the hot case (WB02, WA02) and centrifugal support 
in the colder case (WA05) . Gravity is here the strongest 
force towards the origin and the situation just reflects 
the radial force-equilibrium we have applied for the inlet 
boundary. This is the hydrodynamic regime. 

At the Alfven point, the residual of the pinch- and 
toroidal pressure-force (jp x S^) is the main collima- 
tor, balanced by the centrifugal term. Thermal pres- 
sure quickly looses importance. This is the magneto- 
hydrodynamic regime. 

In the asymptotic region beyond the light-cylinder, de- 
coUimation by electric forces overcomes the centrifugal 
force and is balanced by the poloidal magnetic pressure 
that changes its sign at the light-cylinder (best seen in 
WBOl). This is the relativistic regime. 

To get a global impression on the relative importance 
of the individual forces we show a radial cut through- 
out the asymptotic jet in figure |9] The strongest forces 
arise across the inner asymptotic light surface which sep- 
arates field lines of high angular velocity from those in 
the non-rotating corona along the axis. Here the electric 
de-collimation is essential. The B^{r) profile is curled 
up from the inner disk radius along the outfiow - result- 
ing in a magnetic pressure gradient that works in unison 
with the toroidal field pinch force until at some radius 
the toroidal field surpasses its maximum and decreases 
(negative gradient). 

The strong gradients in toroidal field and rotation in- 
duce a current sheet and give rise to an electric charge. 
The space charge pc — (l/4)7rV • E is positive close to 
the axis and changes its sign at a critical line as defined 
bylGoldreich & Julian|(|1969 1. 



4.2.2. Energy conversion 

In the simulations where injection is sub-magnetoslow, 
the energy flux is not a free parameter, but is consis- 
tently determined by the simulation of the disk-wind. It 
is hence of interest how the partitioning and conversion 
is realized. From the values of /imax given in Tab.|2]it is 
obvious that our disk-corona suppo rts on ly mildly rela- 
tivistic flows below r — 1.5 (section 2.3.21. 

In Fig. [To] (bottom left panel) we show the efficiency 
cr of Poyntmg flux to kinetic flux conversion along the 
field line with rfp = 2 in the fast component of the jet. 
Here, a is below equipartition already at the inlet and 
it further decreases as F approaches ^. In Fig. 10 (top 
left) the toroidal velocity shows that the fiow decouples 
from co-rotation with the magnetic field at the Alfven 



point. Beyond the Alfven point, angular momentum is 
then carried predominantly by the magnetic field. The 
poloidal velocity increases from low injection value (sonic 
velocity) to 0.5c. Further acceleration can not be ex- 
pected as the bulk of the energy is already in kinetic 
form. The right panel of Figure [TO] shows the individual 
energy channels compared to the rest-mass fiux for the 
same field line. At the base of the jet, the strong poloidal 
electric currents (a strong toroidal field) give rise to an 
outfiow with JC <T< —Q < S < M, predominantly 
transporting energy via rest-mass and Poynting-fiux. 

The kinetic energy fiux surpasses the thermal fiux at 
the Alfven point and further overcomes the gravitational 
binding energy term shortly thereafter. This is not 
surprising, since the escape-surface can be close to the 
Alfven surface at least for the inner field lines (see also 
Fig. |7|. 

Only then, the cold limit /i = r(cr + 1) is applicable 
- it is certainly valid in the asymptotical outflow where 
thermal and gravitational energy fluxes are negligible. 

4.3. Dependence on the launching environment 

For the simulations described up to now, we have per- 
formed in addition several parameter runs in order to 
investigate how the resulting jet dynamics depends on 
the (prescribed) launching conditions - the disk corona 
(see Tab. [2|. We now focus on the impact of the plasma 
P and the aisk temperature parameter e. 

In general, a low /3 (a stronger magnetic field) we find 
that the outflow tends to coUimate more, as indicated 
by the higher average coUimation degree ^ and a lower 
momentum weighted jet radius rj^t- This in principle de- 
creases the MHD acceleration efficiency which critically 
depends on the divergence of ffux surfaces. It is straight- 
forward to define the mass fiux-weighted (half-) opening 
angle of outflow. 



6j{, ~ atan ^ 



< 



'M 



(48) 

< 7° for the out- 



which translates to angles of 3° 
flows under consideration. 

The impact of the magnetic fleld strength on the 
amount of mass flux is not clearly visible, as the two 
simulations with e = 1/6, 2/3 show a different trend. As 
the e-parameter is simply a proxy for the disk corona 
density, it will affect the coUimation in the following 
manner: A higher inffow density lowers the Alfven sur- 
face towards the disk surface which in turn broadens the 
current topology and therefore widens the ffow. This is 
also the trend that we observe in the indicators ^ and 
Tjct. Given that the injection speed calculated itera- 
tively from the outffow simulation approaches the slow- 
magnetosonic speed, we expect the mass ffux to scale as 
M oc ^ypp. In fact, this is approximately realized since 

we have M(e = l/6)/M(e = 2/3) = ^ 1-6. 

The change of the initial split-monopole inclination 9 
has little effect on the overall jet coUimation angle. In 
particular, we observe an opposite trend as the wider 
initial field with 9 — 77° ends up slightly more coUimated 
than the one with 9 = 85°. Clearly a wider initial field 
leads to a larger jet radius rjet. Here, we like to stress 
the point that for the final steady state solutions in our 
simulations the initial field structure is important only 
insofar as it also prescribes the poloidal magnetic field 
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Fig. 8. — Accelerating (left column) and coUimating (right column) forces along the field line rooted at Tfp = 2 for the three reference 
simulations (from top to bottom): WBOl (split monopole), WA02 (hourglass potential field, hot case), and WA05 (hourglass potential field, 
cool case). Shown are the contributions from gravity, gas pressure gradient, centrifugal force, poloidal magnetic field pressure gradient, 
poloidal field tension and pressure gradient, toroidal magnetic field pressure gradient and tension, and forces due to the electric field. 
Vertical lines indicate the critical surfaces - the Alfven point along this field line, denoted by 'A', the light cylinder radius denoted by 'Ic', 
and the fast magnctosonic point denoted by 'F'. In this logarithmic representation a change of sign in the force direction is indicated by 
the singularities along the graphs. 



profile along the outflow launching boundary. The field 
structure is completely changed from the initial steady 
structure to a new dynamic equilibrium. Thus it makes 
no sense to compare the coUimation of the initial field 
with the coUimation of the outflow fleld distribution. 

Having pointed out the crucial role of the vertical en- 
ergy flux from the disk surface /i and the closely related 
quantity a = S/{K, + M.), we now study the two most 
promising handles in increasing /i. That is (i) a decrease 



in mass flux M and (ii) an increase in Foynting-flux S. 
We first focus on (i) and describe (ii) thereafter. 

4.3.1. Towards low mass loading 

A way to obtain high-speed jets seems to be a lower 
mass load injected into a similarly-strong m agnetic 
flux. According to the well-known Michel-scaling ( Michel 
19691, the asymptotic outflow velocity depends on the 
mass flux Uaa oc M~^/'^. We investigate this interrelation 
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Fig. 9. — Trans-field force cut at z = 200, inner (Ici) and outer 
(Ico) light-cylinder. The differentially rotating field lines are fastest 
at the inner disk-radius, resulting in the inner light-cylinder - here 
electric dc-coUimation is important. The B^{r) profile is curled 
up from the inner disk radius onwards and results in a magnetic 
pressure gradient that works in unison with the pinch force until 
the toroidal field surpasses its maximum. Within r < 3, we omit 
the curves for thermal and poloidal pressure. These terms fiuctuate 
around ±10"^ while balancing each other, (run WA02) 



running another set of simulations where we change the 
mass flux by prescribing a low injection velocity (instead 
of lowering the density of the injected gas, see Tab. |3| 

However, we observe that for low injection speeds 
Vz < 0.4 W0 the resulting (numerical) mass flux does not 

follow the expected linear relation M cx v-my Instead, for 
lower and lower injection speed it approaches some seem- 
ingly unphysical offset value (Fig. |11[ ). This value is un- 
physical in the sense that it departs from the prescribed 
value at the boundary condition dMiaj/dr = 27rr/5inji'inj- 
This difficulty arises because the injected ffow is sub- 
magnetosonic and thus overdetermined by simultane- 
ously assigning a proflle in p and Vz- 

Although this is in principle problematic, it is not ncc- 
essarely fatal for the investigation of low mass flux out- 
flows. The mass flux is governed by the magneto-slow 
point and is thus re-arranged along the flow. Indeed we 
find that above the magneto-slow surface the field line 
constants are very well conserved. This indicates that 
the ffow dynamics transcends at the magneto-slow sur- 
face from a seemingly unphysical state into a MHD ffow 
that satisffes the critical conditions at the magnetosonic 
surfaces. 

This technique of self-adjusting the sub-slow mass in- 
ffow, however, turned out to be limited towards lower 
mass ffuxcs. The resulting mass ffuxes are unfortunately, 
still too high and do not allow substantially higher mag- 
netisation (e.g. /imax = 2.01 for simulation 1001 com- 
pared to ^max = 1-33 in run WAOl). Figure 11 shows 
the trends concerning coUimation ^, jet radius rjot, and 

integral mass ffux M in the asymptotic ffow. Increasing 
the mass ffux enhances coUimation while /^max decreases 
accordingly. For high injection speed, we observe that 



the wind originating from the very inner disk evolves 
into a thin ballistic ffow layer that has little in common 
with the jets we are interested in. 

4.3.2. Poynting dominated flows 



Given the limitations mentioned above, we prescribe a 
priori the limiting energy ffux parameter /n by constrain- 
ing both the mass ffux and the Poynting ffux along the 
injection boundary - with the hope of thus providing a 
sufficiently energetic disk wind. 

To achieve this, we adopted a fixed-in-time toroidal 
magnetic field distribution, oc —rj/r which necessarily 
changes ri^(r)*. The toroidal field distribution following 
a 1/r profile corresponds to jz = and has a profound 
physical motivation as it anticipates the radial currents 
expected in the disk-corona. 

Table [4] summarizes the simulation runs performed 
within this setup. These simulations have a considerably 
stronger toroidal magnetic field at the injection point. 
We apply up to ~ 4i?p. 

An exemplary process enhancing large-scale toroidal 
fields in the disk corona could be the MRI-driven d ynamo 
under current investigation by many auth ors (e.g. Miller 
fc Stone|[2000| |von Rekowski et al.||2003| . The jet even 



tually evolving trorn these disks is not propelled by the 



Blandford & Payne ( 1982 1 m echanism, but driven by the 
toroidal magne ticpressure f Contopoulos 1996 De Vil 



|Uers et al.|2003||Kato et al.|2004j . 
lets have initially been proposed 



These so called lower 



jets have initially been proposed by Lynden-Bell ( 1996 1 
and directly extract Poynting-fiux from the dist rather 
than first converting rotational energy into the twisted 
magnetosphere that is present at the Alfven surface.^ 

The simulations described here are of a mixed type, 
since they combine large-scale open field lines with a 
toroidal field emerging from a radial current. For an 
example simulation of this type, we show the conversion 
of energy in Fig. [12] similar to Fig. [TOjfor the low-energy 
case. By design, the injected Poynting fiux surpasses the 
rest mass fiux with cr ~ 5 within the fast outfiow com- 
ponent. 

The outflow, which is initially Poynting flux domi- 
nated, does not reach equipartition at the fast magne- 
tosonic s urface r = rp, where we merely find F ~ fj}/^ 
following Michel] ( ]1969] ); ]Beskin et all ( ]1998] ). In the 



asymptotic region r rp, the length scales for addi- 
tional fiow acceleration and coUimation wo uld increase 



expon entially with T oc (/itlnr) ' (e.g. Tomimatsu 
(]1994|), which is clearly beyond the reach of our numer- 
ical method. 

Our simulations indicate that, given sufficient Poynt- 
ing flux, bulk Lorentz factors derived from AGN jet 
observations can be obtained within several hundred 
Schwarzschild radii, F ~ 6 for model M08. Since ac- 
celeration has proven to be most effective around the 
Alfven point which is expected to be very close to the 
central object (rA < He), we expect this conclusion to 
remain valid also for higher /x and thus higher terminal 
F ffows. However, we note that when increasing a, the 
energy conversion efficiency decreases, a situati on com- 
monly denoted as g-problem in puls ar winds (Rees & 



Gunn 1974 Kennel & Coroniti 1984). Several authors 



have recently addressed this issue with partly controver- 

* In case of a central black hole causality requires that rf2^(r) < 
0.6 which corresponds to the ISCO velocity for the Schwarzschild 
case wit h rjqnn = 1 in our scaling. 

^ See Kato for a revi ew. Note also that t hese jets have 

successfully been reproduced by |Lebedev et ar| | |2005[ l in laboratory 
experiments with purely radial current distributions at the base. 
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Fig. 10. — Dynamical quantities as a function of the distance along the field line rooted at rfp = 2 for the model WA02. Vertical lines 
indicate the slow-magnetosonic, Alfvcn, light surface and fast magnetosonic transitions (from left to right). Left: Isorotation parameter 
n^, poloidal and toroidal velocity are shown in the top panel. Lorentz factor F, energy conversion efficiency a, and normalized total energy 
flux fi in the bottom panel. The flow is below equipartition already at the base of the jet. Right: Complete energy flux ratios. In the 
asymptotic region, thermal and gravitational fluxes are obviously negligible. 



TABLE 3 

Effect of mass-loading 



ID 


Top 


/3 


e 




V 


Remarks 


Fmax 






''^p,max 




M 


1001 


A 


0.2 


2/3 


0.01 


var. 




1.66 


2.01 


25.61 


0.75 


21.11 


10.63 


101 


A 


0.2 


2/3 


0.1 


var. 




1.52 


1.76 


25.58 


0.71 


20.49 


13.78 


102 


A 


0.2 


2/3 


0.2 


var. 




1.38 


1.55 


26.02 


0.65 


19.34 


17.53 


104 


A 


0.2 


2/3 


0.4 


var. 




1.31 


1.45 


25.30 


0.60 


17.07 


25.84 


108 


A 


0.2 


2/3 


0.8 


var. 




1.16 


1.24 


40.05 


0.47 


14.74 


51.30 


112 


A 


0.2 


2/3 


1.2 


var. 




1.21 


1.21 


33.77 


0.50 


14.10 


81.88 



Note. — As in table [21 but for the extended parameter study with a given t^inj ■ 



TABLE 4 

POYNTING DOMINATED FLOWS 



ID 


Top 


/9 


e 




V 


Remarks 


Fmax 


/^max 


5 


^p.max 


net 


M 


M02 


A 


0.2 


2/3 


0.1 


2 




2.18 


3.29 


27.55 


0.86 


23.71 


22.48 


M04 


A 


0.2 


2/3 


0.1 


4 




3.51 


7.46 


8.96 


0.94 


29.00 


48.85 


M08 


A 


0.2 


2/3 


0.1 


8 




6.11 


25.61 


6.8377 


0.97 


35.89 


80.40 



Note. — As in table|2j but for the extended parameter study with a given Vi^j and 7]. 



sial results (| Komissarov et al.||2009 Tchekhovskoy et al. 
2009 Lyubarsky 2009 ) , so that atter the successtul accel- 



eration towards relativistic speeds, Poynting-flux could 
still remain and one is tempted to ask: Is there a a- 
problem for AGN-jets? The simulations presented here 
are not fit to answer this question satisfactory. How- 
ever, we certainly know that the mildly relativistic disk 
winds presented earlier do not suffer from this, as they 
are launched already in sub-equipartition. 

5. SUMMARY 

We have presented ideal MHD simulations of the 
formation of special relativistic disk winds using the 
FLUTO 3.0 code. On the technical side, the key points 
are: 

i) The inclusion of (Newtonian) gravity allows us to spec- 
ify an astrophysically sensible boundary condition of a 
hydrodynamically stable disk corona. We can thus con- 
sistently follow the acceleration from initially sub-escape 
velocity winds. 



ii) Much dedication has been put in the development 
and testing of a novel realization for the outflow bound- 
ary that enables us to simulate for hundreds of inner 
disk rotations while minimizing spurious coUimation due 
to artificial boundary currents. Our detailed study of jet 
coUimation is possible only through this effort. 

As a general result we obtain well coUimated jets with 
a mass flux weighted half-opening angle of 3 — 7° and 
mildly relativistic velocities depending on the launching 
conditions for the outflow. The flow coUimation happens 
mainly in the classical (non-relativistic) regime before the 
light surface. A major result of our simulations is that 
we - for the first time - self-consistently calculated the 
shape of that light surface. The light surface determines 
the "relativistic" charater of the flow. Material which 
traverses the light surface experiences the full relativistic 
effects. 

We can identify three dynamically distinct regions in 
terms of flow coUimation. 
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Fig. 11. — Jet coUimation and dynamics against injection speed 
parameter Vi^j, The horizontal line indicates the mass flux when 
Vinj is not specified (WAOl). For iiinj < 0.4, the mass flux ap- 
proaches an obviously unphysical offset value. At higher vi^j , the 
mass flux follows the expected linear dependence on the injection 
parameter (thin solid line). 
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The relatively slow winds found to arrise at large dis- 
tances of around 100 Schwarzschild radii may be ob- 
served as X-ray absorption winds in radio-quiet AGN. 

In the case of Blandford-Payne disk winds (i?^ = ini- 
tially), the outflow is kinetic energy-dominated with the 
ratio of electromagnetic energy flux to kinetic energy flux 
(T < 1 already at the jet base and with this ratio further 
de-creasing downstream the outflow. These disk winds 
start out at sonic speed and reach only mildly relativistic 
speeds up to Lorentz factors F < 1.5. We have also in- 
vestigated cases where the jet Poynting flux is increased 
by directly injecting an additional toroidal magnetic fleld 
from the disk boundary. In this case we achieve terminal 
Lorentz-factors up to Fmax — 6 while the super fast- 
magnetosonic jet remains Poynting-flux dominated. 

We thank Andrea Mignone and the PLUTO team for 
the possibility to use the PLUTO code and for his sup- 
port with the relativistic module. Numerical simula- 
tions were performed on the PIA cluster of the Max 
Planck Institute for Astronomy (Heidelberg) located at 
the Rechen-Zentrum in Garching. O.P. likes to thank 
Bhargav Vaidya for his commitment in numerous discus- 
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Fig. 12. — As in Fig. |10| however calculated for simulation 
run M04. The injected flow is strongly magnetized and remains 
Poynting-dominated when leaving the computational domain. 

i) In the hydrodynamic regime upstream of the Alfven 
surface, gravity balances thermal and magnetic pressure, 
respectively the centrifugal force in the colder case. 

ii) In the magneto- hydrodynamic regime following the 
the Alfven surface downstream, the residuals of mag- 
netic pinch and the toroidal magnetic pressure gradient 
balances the centrifugal force. 

in) In the relativistic regime located downstrem of the 
light surface, the poloidal magnetic pressure gradients 
now impose a coUimating force against electric field de- 
coUimation. Electric forces ultimately overcome the clas- 
sical magneto-centrifugal contribution. 

A steep rotation profile of the field line as given by a 
Keplerian disk results in a light surface geometry which 
steepens for large radii. Depending on the magnetic field 
profile, the light surface may even coUimate along the 
flow for large radii. In such a case the relativistic core 
inside the light surface is naturally confined by a non- 
relativistic wind. The ability of both the relativistic jets 
and the non-relativistic disk winds to coUimate may pro- 
vide confining agents for an axial ultra-relativistic funnel 
which could probably launched by the Blandford-Znajek 
process. 
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APPENDIX 
ZERO CURRENT BOUNDARY 

One of the major goals of this paper is to investigate the colliniation behavior of relativistics outflows. It is therefore 
essential to exclude any numerical artifacts leading to a spurious flow coUimation. We find that the standard zero- 
gradient outflow boundary conditions may lead to an un-physical Lorentz force in radial direction implying such 
spurious coUimation (or de-coUimation) . 

Thus, we put substantial effort in implementing and testing an enhanced outflow boundary condition to the code. 

To get a handle on the Lorentz-force x Bp, one has to address the toroidal electric currents at the grid boundary. 
In principle there are (at least) two options. One is the possibility to copy the toroidal electric current across the 
boundary. While this approach should minimize spurious coUimation efficiently, we observed that the overall stability 
of the simulation was decreased. Thus we decide to use the following zero-toroidal current outflow boundaries in our 
simulations. 

In this case we take advantage of the staggered grid by enforcing zero toroidal currents while simultaneously satisfying 
the solenoidal condition V • B = 0. In the following our procedure is described in detail. We consider computational 
grid cells (iend,j), adjunct to the domain boundary at (icnd + l,j), as illustrated in Fig. 13 The magnetic field 



components of the domain, -Bt(«endii + 1/2), -Bn(*cnd + 1/2, j), Bn{icnd + 1/2, j -I- 1), togetEer with the transverse 
field component -Bt(icnd + 1, j + 1/2) of the first ghost zone, constitute a toroidal corner-centered electric current 
l4>{icnd + 1/2, j + 1/2). Utilizing Stokes theorem, I,p — JdS • V x Bp ~ fdl ■ Bp, we then solve for the unknown field 
component -Bt(*cnd + 1, J -l- 1/2) under the constraint that = 0, 



Ar 



(Al) 



where we have assumed an equally spaced grid for clarity of the argument. Once i3t(iend + I7J + 1/2) is known for 
all J, the next layer of normal field components i?n(*cnd + 3/2, j) can be inferred from the V • B = constraint in its 
integral form. 



5nk„i+3/2j + l — 



3+3/2) 



A5„ 



io„d+3/2J+l 



(A2) 



For the next grid layer, the transverse field components can again be found applying Eq. |A1| and the process is 
repeated for the each layer. 
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Fig. 13. — Construction of the V X Bp = and V ■ B = boundary condition. Shown is the last grid slab of the domain (iendii) ^^nd a- 
ghost zone of two elements. 



Some words of caution. We find that the current-free magnetic field boundary condition can only be realized when 
the grid cell aspect ratio Az/Ar is not too large. An aspect ratio of e.g. 12/1 resulted in errors of 100% in B^ at the 
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most critical areas close to the symmetry axis leading to an overall unstable flow evolution. We flnd that as a rule of 
thumb, an aspect ratio of 3/1 should not be exceeded. We also emphasize that it is essential to treat the grid corners 
consistently. This is because field components in the corner, i3t(iend + 1, jcnd + 1/2), and -Bn(*cnd + 1/2, jond + 1), are 
interrelated which would lead to an ambiguity. In order to avoid this ambiguity, we decided to extrapolate the values 
in question which docs provide the information that is missing otherwise. 

We demonstrate quality of our approach by showing results of simulations which do not apply the zero current 
but the zero-gradient or the zero second derivative outflow condition with otherwise the same flow parameters as in 
simulation WA04 (Fig. 14 1. As it can be seen, for zero-gradient boundary conditions, the effect of coUimation by 



artificial currents is so strong that no steady-state can be reached and the flow is continously squeezed towards the 
axis. Also in zero second derivative, we observe an artificial alignment with the grid geometry. 
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Fig. 14. — Comparison of simulations applying a variation of outflow boundary conditions for the magnetic fields at the time of 100 inner 
disk rotations. The parameters are equal to those in simulation WA04. The grayscale indicates the Lorentz factor log(r — 1) as in figure 
W] the poloidal magnetic field (the poloidal electric current) is shown in thick (thin) white contours. Standard zero gradient, zero second 
oerivative, zero current boundary condition, respectively (from left to right). 



We check the geometry dependence of the zero current outflow boundary by several r ealiz ations of the fiducial 
run WA04, each with the same resolution but with a different grid size or shape. Figure [Ts] (left panel) compares 
the steady state flow characteristics for various boxes with ratios Az/Ar G {1/1,2/1,4/1}. While geometries and 
sizes with Az/Ar > 2/1 are in excellent agreement, the quadratic domains show significantly thinner characteristics. 
The reason for this discrepancy is the sub Alfvenic flow that traverses the Z^^d boundary in "broad" domains. In 
these underdetermined simulations, current circuits start to unclose at the sub Alfvenic part of Zend which ultimately 
destro ys the Butterfly shape in the entire domain. As also noticed and extensively discussed by |Krasnopolsky et aL] 



(19991, a sub Alfvenic (vertical) outflow can not obtain the proper critical point information and leads to erroneous 
extensive coUimation. This problem can be avoided by taking the position of the critical Alfven surface into account, 
hence we choose a ratio of 2/1 for our science simulations. 
Finally we check convergence by comparison to a half-resolution run with 256 x 5 12 grid elements. The solutions are 



in good agreement, indicated by contours of the Alfven mach number in figure 15 (right panel). In conclusion we use 
a grid of 512 x 1024 cells with a domain size of (r, z) = (102, 204) inner disk radn7 ensuring that the presented results 
depend mostly on the disk corona boundary. 

THE INJECTION BOUNDARY IN MHD-JET SIMULATIONS 

The number of constraints imposed on a giv en boundary mus t equal the number of waves allowed to travel through 
the boundary to the simulation domain (e.g. 'Bogovalov' (1997)). In MHD the number of characteristics equals the 
number of "variables" minus one - due to the V • B = constraint - to give a total number of seven. 

Another way of looking at this is by considering the critical surfaces as internal boundaries. For example, in a 
sub-sonic flow, the intersection of the two characteristics C± with respective wave velocities w ± Cj determines the 
hydrodynamical state. With respect to any sub-sonic boundary condition, only one characteristics is incoming, while 
the outgoing characteristics originates at the sonic point where its velocity vanishes. The same is true for the supersonic 
case, only now C_ transports the information of the sonic point downstream. 

Applying a, x — t diagram one may understand why the sonic point constitutes a fixed-in-time boundary condition. 
That is because here C_ becomes singular, and hence this part of information ( the Riemann-invariant) starts to travel 



up- and down-stream from the critical point (see also Landau & Lifshitz 



Following these general considerations, we see that the correct number o: 



1959 chap. X). 



constraints for a sub-magnetoslow boundary 



is four, since the flow is expected to pass through three characteristics. In other words, there are four outgoing waves: 
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Fig. 15. — Comparison of different grid realisations with zero current outflow boundaries after 200 inner disk rotations. Parameters as 
in simulation WA04. Left: Various geometries and sizes. Shown is the critical Alfven surface and the light cylinder. Right: Convergence 
test with two grid resolutions. Shown are contours of the Alfven mach number M for 256 X 512 and for 512 X 1024 grid elements. 

the slow magnetosonic wave, the entropy wave, the Alfven wave and the fast magnetosonic wave. Naturally, along a 
boundary where the flow is super magneto-slow this number equals five. Within this limited freedom, those boundary 
conditions which best prescribe the astrophysical problem should be used. 

To allow the outflow to settle i nto a steady state requires certain conditions to be met also at the boundary. In our 
paper, we follow the argument by Krasnopolsky et al. ( 1999 1. According to the (axially symmetric) induction equation 
it is 

dtB,^l/rdrE^; dtBr ^ -d,E^- dtB^ = d,Er - drE,. (Bl) 

Since in steady state dtB^ = dtBr = 0, the only physical solution is E^f, = which is satisfied by Vp||Bp. 

A number of authors prescribe in addition a fixed-in-time value for Er = ^l^ Bz, however, this is equivalent to keeping 
the iso-rotation parameter f2^ constant in time, as the time evolution of B^ is already suppressed by the choice of E^. 
Unlike often stated, the certainly proper choice of constraining {E^,Er) is not dictated by the ideal MHD condition - 
vanishing electric fields in the co-moving frame - (which is respected by design) , but by the steady-state considerations 
given above. In the literature of "disk-as-boundary" jet formation simulations a variety of choices for the injection 
boundary exist. Table |5] reviews a couple of them in chronological order. 
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TABLE 5 

Injection conditions for. disk-as boundary simulations 



Authors 



sub/super slow 



Outgoing waves 



Constraints 



Nature of constraints 



Ustyugova et al.l ||1995|l'' 


sub 


4 


4 




Ouycd Pudritz ('19971'' 


sub/super 


4 


6 


E^,v^, Vp, B^, p, s 


Komanoya ct ai. ( 1997T* 


sub 


3 


3 


E^,Er,Vz 


ITrasnopolsltv ct al. (I999p 


super 


5 


5 


E^,Er,p,S,Vz 


Ustyugova ct al. ( i999|^ 


sub 


4 - 5 


4 


E^,Er,p,S 


Komissaroy et al.fKODVp 


super 
sub - super 


5 

4 - 5 


5 

4 - 6 


E^,B'^,p,v'i,n'' 
E^,n^,p,s 


'I'his worit'* 



Ouyed fc Pudritz| ( |1997 i TlFendt fc Cemeljic (20021 allow no feedback from the jet to the disk and seem to over-determine 
;ir simulations. We were able to reproduce the consequences with our own simulations. A numerical boundary layer develops 



b 

their simulations. We were able to reproduce the consequences 
which creates a steep gradient 9zp(r, z) as the code tries to match the dense disk-boundary with the jet-solution (their figure 
4). To conserve mass flux, the poloidal velocity will jump within a few grid cells, a spurious acceleration which is independent 
of resolution. Due to their additional Alfven pressure, the injection velocity is not entirely clear to us, but we suspect it to be 
d ynamically sub-slo w. 



Romanova et al 



1997) perform isothermal simulations without solvin g the energy equation. Since the sub-slow injection 



constraints on the following three quantities are given, v^, E^, Er - just as Ustyugova et al. ( 1995 1. 



inj ect ion speed is a" 
f jKomissarov et al. (20071 also inject super-slow 



The 



Krasnopolsky ct al. ( 1999) prescribe mass flux with the choice {E^, Er, s, p, Vz) and super-slo w injection. 

Ustyugova et al.^ (,1999 J prescribe density p instead of the vertical velocity Vz (thus change from Ustyugova et al. ( 1995 1). 
Jlo wed to become super-sonic which strictly speaking results in an underdetermined system. 

As the energy equation is not solved, one condition is already used for 
fixing the entropy. Constraints in the injection boundary *far from the disk) disk are (-E^, B'^ , p, w'', fl^) where the 7;-coordinate 
describes a "radial" direction (r^ 
s See the discussion in section [37 



— /a + z'^ , a > 1) in elliptical coordinates. 
Lljof this paper. 
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